script by: Michael T. Hallworth

#############################################################################
#
#   System set up (uncomment code if running code for the first time)
#
#############################################################################

# Alternatively you can run the following few lines of code to get the packages
# library(devtools)
# install_github("SLisovski/GeoLight", ref = "Update_2.01", force = T)
# install_github("SLisovski/GeoLocTools")
# library(GeoLocTools)
# setupGeolocation()
# devtools::install_github("MTHallworth/SGAT")

Load required packages

################################################################################
#
#                           Load packages 
#
################################################################################

library(GeoLocTools)
data(wrld_simpl, package = "maptools") 

setupGeolocation()

library(SGAT)
library(raster)
library(sp)
library(landscapemetrics)
library(LLmig)
library(rcartocolor)

Read in spatial layers

# Read in spatial layers 
Americas <- raster::shapefile("Spatial_Layers/Americas.shp")

OSFLDist <- raster::shapefile("Spatial_Layers/OliveSidedFlycatcher.shp")

States <- raster::shapefile("Spatial_Layers/st99_d00.shp")

Canada <- raster::shapefile("Spatial_Layers/Canada.shp")

## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum D_North_American_1983 in Proj4 definition: +proj=aea
## +lat_0=40 +lon_0=-96 +lat_1=50 +lat_2=70 +x_0=0 +y_0=0 +ellps=GRS80 +units=m
## +no_defs

AK <- subset(States, NAME == "Alaska")

#############################################################################
#
#                          Define projections used in the analysis
#
#############################################################################

WGS84 <- sf::st_crs(4326)$proj4string
NAEA <- "+proj=aea 
         +lat_1=20 
         +lat_2=60 
         +lat_0=40 
         +lon_0=-96 
         +x_0=0 +y_0=0 
         +ellps=GRS80 +datum=NAD83 +units=m +no_defs"

# Transform Canada into WGS84

Canada <- spTransform(Canada,CRS(WGS84))

# Set the projection for shapefiles
crs(States) <- WGS84
crs(Americas) <- WGS84
crs(OSFLDist) <- WGS84

################################################################################
#
#                          Read in light-level data 
#
################################################################################

AncFiles <- list.files(path = "Data/Anchorage_recoveries",
                       pattern = ".lux",
                       full.names = TRUE)

FairFiles <- list.files(path = "Data/Fairbanks_recoveries",
                        pattern = "driftadj.lux",
                        full.names=TRUE)

TetlinFiles <- list.files(path = "Data/Tetlin_recoveries",
                          pattern = ".lux", 
                          full.names = TRUE)

# Read just the file names for Bird ID

AncNames <- list.files(path = "Data/Anchorage_recoveries",
                       pattern=".lux")

FairNames <- list.files(path = "Data/Fairbanks_recoveries",
                        pattern="driftadj.lux")

TetlinNames <- list.files(path = "Data/Tetlin_recoveries",
                          pattern = ".lux")

# Combine the files into a single object

OSFLFiles <- c(AncFiles, FairFiles, TetlinFiles)

# Combine the bird ID
BirdId <- c(AncNames, FairNames, TetlinNames)

# More readable BirdId 
BirdId <- substr(x = BirdId,
                 start = 1,
                 stop = 4)

#Determine the number of birds
nBirds <- length(BirdId)

## ----readLux-------------------------------------------------------------
OSFLdata <- lapply(X = OSFLFiles, 
                   FUN = readMTlux) 

for(i in 1:nBirds){
  #Take the log of light data#
  OSFLdata[[i]]$Light <- log(OSFLdata[[i]]$Light)
}

# Take a quick peek #
# str(OSFLdata,2)

###########################################################################
#
#                          Read in Capture Locations 
#
###########################################################################

## ----caplocs-------------------------------------------------------------
# Set the capture coordinates for each bird #
CapLocs<-array(NA,c(nBirds+1,2))

#CapLocs[1,] <- c(-149.60425,61.38084)    # F202
#CapLocs[2,] <- c(-149.74919,61.28071)    # F203
#CapLocs[3,] <- c(-149.80309,61.29737)    # F311
#CapLocs[4,] <- c(-146.832960,65.335090)  # K288 
#CapLocs[5,] <- c(-147.1002,64.71122)     # K648
#CapLocs[6,] <- c(-146.800030,65.33139)   # K654
#CapLocs[7,] <- c(-146.943180,65.34078)   # K657
#CapLocs[8,] <- c(-142.30559,63.13592)    # K659
#CapLocs[9,] <- c(-146.832960,65.335090)  # K660
#CapLocs[10,] <- c(-147.1002,64.71122)    # K665
#CapLocs[11,] <- c(-147.758300,64.936900) # K670
#CapLocs[12,] <- c(-146.832960,65.335090) # S050 
#CapLocs[13,] <- c(-147.00492,64.78862)   # S052 
#CapLocs[14,] <- c(-142.30217,63.13515)   # Q330 
#CapLocs[15,] <- c(-143.21355,63.16566)   # Q334
#CapLocs[16,] <- c(-148.94470,65.37691) # Q335
#CapLocs[17,] <- c(-146.832960,65.335090)  # K288 

CapLocs[1,] <- c(-149.603083,61.380549)    # F202
CapLocs[2,] <- c(-149.74919,61.28071)    # F203
CapLocs[3,] <- c(-149.80309,61.29737)    # F311
CapLocs[4,] <- c(-148.94351,65.37703)  # K288 
CapLocs[5,] <- c(-148.90802,65.3997)     # K648
CapLocs[6,] <- c(-146.80003,65.33139)   # K654
CapLocs[7,] <- c(-146.943180,65.34078)   # K657
CapLocs[8,] <- c(-142.30559,63.13592)    # K659
CapLocs[9,] <- c(-146.83296,65.335090)  # K660
CapLocs[10,] <- c(-147.1002,64.71122)    # K665
CapLocs[11,] <- c(-147.758300,64.936900) # K670
CapLocs[12,] <- c(-146.8383,65.33454) # S050 
CapLocs[13,] <- c(-147.00492,64.78862)   # S052 
CapLocs[14,] <- c(-142.30217,63.13515)   # Q330 
CapLocs[15,] <- c(-143.21355,63.16566)   # Q334
CapLocs[16,] <- c(-148.94470,65.37691) # Q335
CapLocs[17,] <- c(-148.9447,65.37691)  # K288 

CapLocs <- round(CapLocs, digits = 2)

#############################################################################
#
#                          First glimpse at the light-data 
#
#############################################################################

# ## ----echo = FALSE--------------------------------------------------------
for (i in 1:nBirds) {
  if (nrow(OSFLdata[[i]]) > 1) {
    lightImage(
      OSFLdata[[i]],
      offset = 19,
      zlim = c(0, 10),
      main = BirdId[i],
      dt = 300
    )
    
    # Add solar lines for capture location
    TwGeos::tsimageDeploymentLines(
      date = seq(
        min(OSFLdata[[i]]$Date, na.rm = TRUE),
        
        max(OSFLdata[[i]]$Date, na.rm = TRUE),
        "day"
      ),
      zenith = 95,
      lon = CapLocs[i, 1],
      lat = CapLocs[i, 2],
      offset = 19,
      lwd = 3,
      col = rgb(0, 0, 255, 150, max = 255)
    )
    Sys.sleep(2)
  } else{
    next
  }
}

# onBird <- data.frame(Deploy = rep(as.POSIXct("2011-01-01"),nBirds),
#                     Recover = rep(as.POSIXct("2011-01-01"),nBirds))

#for(i in 1:nBirds){
#lightImage(OSFLdata[[i]], 
#           offset = 19, 
#           zlim = c(0,10),
#           main = BirdId[i],
#           dt = 300)
#onBird[i,1] <- as.POSIXct(locator(n=1)$x,origin = "1970-01-01")
#onBird[i,2] <- as.POSIXct(locator(n=1)$x,origin = "1970-01-01")
#}

# saveRDS(onBird,"onBird_dates.rds")

onBird <- readRDS("onBird_dates.rds")

tempdata <- vector('list', nBirds + 1)
# SUBSET OUT THE DATA TO ONLY INCLUDE WHEN TAGS WERE ON BIRDS #
for (i in 1:nBirds) {
  tempdata[[i]] <-
    OSFLdata[[i]][OSFLdata[[i]]$Date > onBird[i, 1] &
                    OSFLdata[[i]]$Date < onBird[i, 2], ]
  
  # K288 has two years of data in the same file - 
  # here we separate them into '2' birds
  if (BirdId[i] == "K288") {
    tempdata[[nBirds + 1]] <-
      OSFLdata[[i]][OSFLdata[[i]]$Date > onBird[i, 2], ]
  }
}

OSFLdata <- tempdata

###############################################################################
#
#                          Assigning twilight events
#
###############################################################################

# This new way of defining the seed basically creates a vector of all the days
# and adds 8:00 to it that way every dark phase is captured.

seed <- lapply(
  OSFLdata,
  FUN = function(x) {
    seed <- seq(
      from = as.POSIXlt(
        paste(strptime(x[1, 1],
                       format = "%Y-%m-%d", tz = "GMT"), "08:00:00"),
        format = "%Y-%m-%d %H:%M:%S",
        tz = "GMT"
      ),
      to = x[nrow(x), 1],
      by = "day"
    )
    return(seed)
  }
)

# use the 'known' dark periods to assign the twilights #

twl <- mapply(
  x = OSFLdata,
  y = seed,
  FUN = function(x, y) {
    findTwilights(
      tagdata = x,
      threshold = 1.5,
      include = y,
      dark.min = 10
    )
  },
  # 0.75 hours minimum dark period
  SIMPLIFY = FALSE
) # return a list

# take a quick look at the light image with the assigned twlights added on #

for (i in 1:(nBirds + 1)) {
  lightImage(
    OSFLdata[[i]],
    offset = 19,
    zlim = c(0, 10),
    main = BirdId[i],
    dt = 300
  )
  
  # Add assigned twilight times to the light image #
  TwGeos::tsimagePoints(
    date = twl[[i]]$Twilight,
    offset = 19,
    pch = 19,
    cex = 0.8,
    col = ifelse(
      twl[[i]]$Rise,
      rgb(255, 0, 0, 150, max = 255),
      rgb(0, 0, 255, 150, max = 255)
    )
  )
  Sys.sleep(2)
}

##########################################################################
#
#                          Edit twilights
#
##########################################################################
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  # # # #
#
# Edits using preprocesslight
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  # # # # #

# twlProc <- vector('list',nBirds+1)

#for(i in 1:(nBirds+1)){

#cat("\nNEW BIRD! \nprocessing bird ", BirdId[i])
#twlProc[[i]] <- preprocessLight(tagdata = OSFLdata[[i]],
#                          threshold = 1.5,
#                           lmax = 5,
#                           offset = 19,
#                            twilights = twl[[i]])
#}

# saveRDS(twlProc,"ProcessedTwilights_Threshold1.5_Feb2021.rds")

twlProc2 <- readRDS("ProcessedTwilights_Threshold1.5_Feb2021.rds")

# remove edited twlights for analysis #
twlProc3 <- lapply(twlProc2,
                   function(x) {
                     if ("Deleted" %in% colnames(x)) {
                       z <- x[!x$Deleted, ]
                     } else{
                       z <- x
                     }
                     return(z)
                   })

###############################################################################
#
# Some Zeniths are a bit unrealistic - 89 for example.
# For those we'll use the calibration data from the anchorage birds
#
##############################################################################
# Create a vector with the dates known to be at deployment #
calibration.dates <- vector('list', 3) # 3 birds are from anchorage

for (i in 1:3) {
  calibration.dates[[i]] <- c(OSFLdata[[i]][1, 1],
                              as.POSIXct("2013-07-30", tz = "GMT"))
}

## ------------------------------------------------------------------------
# Extract twilight data during calibration period
calibration.data <- vector('list', 3)

for (i in 1:3) {
  calibration.data[[i]] <- subset(
    twlProc3[[i]],
    twlProc3[[i]]$Twilight >= calibration.dates[[i]][1] &
      twlProc3[[i]]$Twilight <= calibration.dates[[i]][2]
  )
}

## ------------------------------------------------------------------------
# create empty vectors to store data #
sun <-
  z <-
  zenith0 <-
  zenith1 <- twl_t <- twl_deviation <- alpha <- fitml <-
  vector("list", 3)

# loop through each of the two individuals #
for (i in 1:3) {
  # Calculate solar time from calibration data
  sun[[i]]  <- solar(calibration.data[[i]][, 1])
  
  # Adjust the solar zenith angle for atmospheric refraction
  z[[i]] <- refracted(zenith(
    sun = sun[[i]],
    lon = CapLocs[i, 1],
    lat = CapLocs[i, 2]
  ))
  
  twl_t[[i]]   <- twilight(
    tm = calibration.data[[i]][, 1],
    lon = CapLocs[i, 1],
    lat = CapLocs[i, 2],
    rise = calibration.data[[i]][, 2],
    zenith = quantile(z[[i]], probs = 0.5)
  )
  
  # Determine the difference in minutes from when the
  # sun rose and the geolocator said it rose
  twl_deviation[[i]] <- ifelse(calibration.data[[i]]$Rise,
                               as.numeric(difftime(calibration.data[[i]][, 1],
                                                   twl_t[[i]],
                                                   units = "mins")),
                               as.numeric(difftime(twl_t[[i]],
                                                   calibration.data[[i]][, 1],
                                                   units = "mins")))
  
  twl_deviation[[i]] <-
    subset(twl_deviation[[i]], twl_deviation[[i]] >= 0)
  
  # Describe the distribution of the error
  fitml[[i]] <- fitdistr(twl_deviation[[i]], "log-Normal")
  # save the Twilight model parameters
  alpha[[i]] <- c(fitml[[i]]$estimate[1], fitml[[i]]$estimate[2])
}

## ----echo=FALSE----------------------------------------------------------
meanAlpha1 <- mean(c(alpha[[1]][1], alpha[[2]][1], alpha[[3]][1]))
meanAlpha2 <- mean(c(alpha[[1]][2], alpha[[2]][2], alpha[[3]][2]))

ALPHA <- alpha[[1]]
ALPHA[1] <- meanAlpha1
ALPHA[2] <- meanAlpha2

# use mean Anchorage sun angle for Fairbanks birds
AncZ <- quantile(c(z[[1]], z[[2]], z[[3]]), probs = 0.95)

## Land is invalid

##############################################################################
#
#                          Use the grouped model to estimate locations
#
##############################################################################
# Convert the light data into something GeoLight can read #
twl.gl <- lapply(twlProc3, export2GeoLight)

# Use the twilight times to distinguish between potential locations #
# you may want to fiddle with the quantile a little bit. The higher the value #
# the fewer the locations #

CLight <-
  lapply(
    twl.gl,
    changeLight,
    quantile = 0.92,
    days = 2,
    summary = FALSE,
    plot = T
  )

# the following edits needed to be made for the mergeSites 
# function to work properly without error
CLight[[7]] <-
  changeLight(
    twl.gl[[7]],
    quantile = 0.90,
    days = 2,
    summary = FALSE,
    plot = TRUE
  )

CLight[[11]] <-
  changeLight(
    twl.gl[[11]],
    quantile = 0.90,
    days = 2,
    summary = FALSE,
    plot = TRUE
  )

CLight[[15]] <-
  changeLight(
    twl.gl[[15]],
    quantile = 0.92,
    days = 2,
    summary = FALSE,
    plot = TRUE
  )

ls_merge <- vector('list', (nBirds + 1))

# # slightly altered GeoLight mergeSites2 function to 
# allow for higher latitudes than 75 N.
#  source('GL_mergeSites_changed.R')
# (a <- Sys.time())
# for(i in 1:(nBirds+1)){
#  cat("Starting analysis for bird # ",i," ",BirdId[i],"\n")
#  ls_merge[[i]] <- mergeSites3(tFirst = twl.gl[[i]]$tFirst,
#                              tSecond = twl.gl[[i]]$tSecond,
#             type = twl.gl[[i]]$type,
#             site = CLight[[i]]$site,
#                         distThreshold = 250,
#             degElevation = -5.28, # AncZ-90 = -5.28
#             alpha = c(2.6,0.9),
#             plot = FALSE,
#             mask = "land",
#             method = "log-norm",
#             spatial_mask = Land)
# print(Sys.time()-a)
# saveRDS(ls_merge[[i]], 
# paste0("D:/OSFL_Geolocator/Merged_sites_",BirdId[i],"_OSFL_Mar2021.rds"))
# }
# Sys.time()

ls_files <-
  list.files(
    "D:/OSFL_Geolocator",
    pattern = glob2rx("Merged_sites*_OSFL_Mar2021.rds"),
    full.names = TRUE
  )

# re-order the saved files so they are in the same order as the other data #
ls_files <- ls_files[c(1:4, 6:12, 16:17, 13:15, 5)]

ls_merge <- lapply(ls_files, readRDS)

twlProc4 <- grouped_twl <- vector('list', (nBirds + 1))

for (i in 1:(nBirds + 1)) {
  Twilight_back <- as.POSIXct(c(ls_merge[[i]]$twl$tFirst,
                                ls_merge[[i]]$twl$tSecond))
  
  Rise <- c(
    ifelse(ls_merge[[i]]$twl$type == 1, TRUE, FALSE),
    ifelse(ls_merge[[i]]$twl$type == 1, FALSE, TRUE)
  )
  
  Site <- ls_merge[[i]]$site
  
  twlProc4[[i]] <- data.frame(Twilight = Twilight_back,
                              Rise     = Rise,
                              Site     = rep(Site, 2))
  
  twlProc4[[i]] <- twlProc4[[i]][!duplicated(twlProc4[[i]]$Twilight), ]
  twlProc4[[i]] <- twlProc4[[i]][order(twlProc4[[i]]$Twilight), ]
  
  grouped_twl[[i]] <- rep(FALSE, nrow(twlProc4[[i]]))
  grouped_twl[[i]][twlProc4[[i]]$Site > 0] <- TRUE
  
  initalGroups <- makeGroups(grouped_twl[[i]])
  
  # Add groups back to twl file
  twlProc4[[i]]$group <- initalGroups
}

behavior <- stationary <- sitenum <- vector('list', (nBirds + 1))
for (i in 1:(nBirds + 1)) {
  # Add behavior vector
  behave <- c()
  for (j in 1:max(twlProc4[[i]]$group)) {
    behave <- c(behave, which(twlProc4[[i]]$group == j)[1])
  }
  
  stationary[[i]] <- grouped_twl[[i]][behave]
  sitenum[[i]] <- cumsum(stationary[[i]] == T)
  sitenum[[i]][stationary[[i]] == F] <- 0
}

############################################################################

## ----setTols-------------------------------------------------------------
tolvalues <- array(NA, c((nBirds + 1), 2))
tolvalues[, 1] <- 0
tolvalues[, 2] <- 0.08

# Manual adjustments for Fall Equinox period
tolvalues[1, 1] <- 0.13
tolvalues[2, 1] <- 0.101
tolvalues[3, 1] <- 0.2
tolvalues[4, 1] <- 0.17
tolvalues[5, 1] <- 0.22
tolvalues[6, 1] <- 0.2
tolvalues[7, 1] <- 0.168
tolvalues[8, 1] <- 0.195
tolvalues[9, 1] <- 0.105
tolvalues[10, 1] <- 0.2
tolvalues[11, 1] <- 0.115
tolvalues[12, 1] <- 0.13
tolvalues[13, 1] <- 0.165
tolvalues[14, 1] <- 0.185
tolvalues[15, 1] <- 0.215
tolvalues[16, 1] <- 0.2
tolvalues[17, 1] <- 0.23

# Manual adjustments for Spring Equinox period
tolvalues[1, 2] <- 0.275
tolvalues[2, 2] <- 0.239
tolvalues[3, 2] <- 0.2
tolvalues[4, 2] <- 0.23
tolvalues[5, 2] <- 0.22
tolvalues[6, 2] <- 0.27
tolvalues[7, 2] <- 0.2285
tolvalues[8, 2] <- 0.2
tolvalues[9, 2] <- 0.2
tolvalues[10, 2] <- 0.162
tolvalues[11, 2] <- 0.188
tolvalues[12, 2] <- 0.17
tolvalues[13, 2] <- 0.22
tolvalues[14, 2] <- 0.24
tolvalues[15, 2] <- 0.2
tolvalues[16, 2] <- 0.2
tolvalues[17, 2] <- 0.22

# devtools::install_github("MTHallworth/SGAT")
library(SGAT)

# Inital path
path <- suppressWarnings(mapply(
  x = twlProc4,
  y = tolvalues,
  FUN = function(x, y) {
    thresholdPath(
      twilight = x$Twilight,
      rise = x$Rise,
      #zenith = x$zenith,
      zenith = AncZ,
      tol = y
    )
  },
  SIMPLIFY = FALSE
))

for(i in 1:(nBirds + 1)) {
  plot(path[[i]]$x[, 2] ~ path[[i]]$x[, 1], main = BirdId[i], pch = 19)
  plot(Americas,
       add = TRUE,
       col = "gray60",
       border = "white")
  points(path[[i]]$x[, 2] ~ path[[i]]$x[, 1], pch = 19)
  points(path[[i]]$x[, 2] ~ path[[i]]$x[, 1],
         type = "l",
         lty = 2,
         lwd = 0.25)
}

# Fixed Locations - I'm going to set the first 10 days as fixed -
# this can be altered #

x0 <- lapply(path, function(x) {
  x$x
})

x0_group <- vector('list', (nBirds + 1))

for (i in 1:(nBirds + 1)) {
  x0_group[[i]] <-
    cbind(
      tapply(path[[i]]$x[, 1], 
            INDEX = twlProc4[[i]]$group,
            median, na.rm = TRUE),
      tapply(path[[i]]$x[, 2], 
             INDEX = twlProc4[[i]]$group,
             median, na.rm = TRUE)
    )
}

z0 <- lapply(x0, trackMidpts)
z0_group <- lapply(x0_group, trackMidpts)

##############################################################################
#
#
#                                Generating a mask 
#
#
#############################################################################
## Function to construct a land/sea mask
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
# set xlim and ylim values need to span the range of your dataset
xlim <- c(-170,-60)
ylim <- c(-89, 90)

abundance <- raster::raster("Spatial_Layers/OSFL_abundance.tif")
binaryOSFL <- reclassify(abundance, matrix(
  c(-1, 0.0001, 0,
    0.0001, 20, 1),
  byrow = TRUE,
  nrow = 2,
  ncol = 3
))

OliveDist <- rasterToPolygons(binaryOSFL, dissolve = TRUE)
OliveDist <- OliveDist[2, ]

STEMprob <- abundance#/cellStats(abundance,sum,na.rm = TRUE)

STEMpts <- rasterToPoints(STEMprob)

# STEMpts[STEMpts[,3]>0,3] <- 2

## ------------------------------------------------------------------------
## Function to construct a land/sea mask
distribution.mask <- function(xlim,
                              ylim,
                              res = c(0.25, 0.25),
                              land = TRUE,
                              shape,
                              values = 1) {
  r <- raster(
    res = res,
    xmn = xlim[1],
    xmx = xlim[2],
    ymn = ylim[1],
    ymx = ylim[2],
    crs = proj4string(shape)
  )
  r <-
    cover(
      rasterize(
        shape,
        shift = c(-360, 0),
        r,
        values,
        silent = TRUE
      ),
      rasterize(shape, r, values, silent = TRUE),
      rasterize(elide(shape,
                      shift = c(360, 0)), r, values, silent = TRUE)
    )
  #r <- as.matrix(is.na(r))[nrow(r):1, ]
  r <- as.matrix(r)[nrow(r):1, ]
  #if (land)
  # r <- !r
  xbin <- seq(xlim[1], xlim[2], length = ncol(r) + 1)
  ybin <- seq(ylim[1], ylim[2], length = nrow(r) + 1)
  
  function(p) {
    r[cbind(.bincode(p[, 2], ybin), .bincode(p[, 1], xbin))]
  }
}

## ------------------------------------------------------------------------
## Define mask for OSFL
STEMmask <-
  distribution.mask(
    shape = SpatialPoints(cbind(STEMpts[, 1:2])),
    values = STEMpts[, 3],
    xlim = xlim,
    ylim = ylim,
    res = c(0.25, 0.25)
  )


## ------------------------------------------------------------------------
log.prior <- function(p) {
  f <- STEMmask(p)
  ifelse(is.na(f) | f == 0, -1000, f)
}

###############################################################################
###############################################################################
#
#
#                                 Fit the inital  model
#
#
##############################################################################
##############################################################################

FIT <- vector('list', (nBirds + 1))

names(FIT) <- c(BirdId, paste0(BirdId[4], "yr2"))

data(wrld_simpl)

# Here you can set the number of iterations to run - 5000 or so
# and the number of samples to thin
niters = 1000
nthin = 5

for (i in 1:(nBirds + 1)) {
  cat("Creating mask for bird ", BirdId[i], ".........\n")

  
  cat("Starting analysis for bird # ", i, " ", BirdId[i], "\n")
  #model <- groupedThresholdModel(twilight = twlProc4[[i]]$Twilight,
  model <- thresholdModel(
    twilight = twlProc4[[i]]$Twilight,
    rise = twlProc4[[i]]$Rise,
    #group = twlProc4[[i]]$group,
    twilight.model = "ModifiedLogNormal",
    alpha = ALPHA,
    # beta = c(2.0, 0.1),
    beta = c(0.7, 0.02),
    x0 = x0[[i]],
    # path
    #x0 = x0_group[[i]],
    z0 = z0[[i]],
    # middle points between the x0 points
    #z0 = z0_group[[i]],
    zenith = AncZ,
    logp.x = log.prior
  )
  
  # define the error shape
  x.proposal <-
    mvnorm(S = diag(c(0.0025, 0.0025)), n = nlocation(model$x0))
  z.proposal <-
    mvnorm(S = diag(c(0.0025, 0.0025)), n = nlocation(model$z0))
  
  cat("Running first model for ", BirdId[i], ".........\n")
  
  # Fit the model
  fit <-
    SGAT::estelleMetropolis(
      model,
      x.proposal,
      z.proposal,
      iters = niters,
      thin = nthin,
      chains = 3
    )
  
  zsum <- locationSummary(fit$z)
  xsum <- locationSummary(fit$x)
  
  proposal.x <- mvnorm(S = diag(c(0.0025, 0.0025)),
                       n = nlocation(cbind(xsum$'Lon.50%', xsum$'Lat.50%')))
  proposal.z <- mvnorm(S = diag(c(0.0025, 0.0025)),
                       n = nlocation(cbind(zsum$'Lon.50%', zsum$'Lat.50%')))
  
  cat("Fine tuning the model for ", BirdId[i], ".........\n")
  
  fit <- estelleMetropolis(
    model = model,
    proposal.x = proposal.x,
    proposal.z = proposal.z,
    x0 = cbind(xsum$'Lon.50%', xsum$'Lat.50%'),
    z0 = cbind(zsum$'Lon.50%', zsum$'Lat.50%'),
    iters = niters,
    # This value sets the number of iterations to run
    thin = nthin,
    chains = 3
  )
  
  
  zsum <- locationSummary(fit$z)
  xsum <- locationSummary(fit$x)
  
  proposal.x <- mvnorm(S = diag(c(0.0025, 0.0025)),
                       n = nlocation(cbind(xsum$'Lon.50%', xsum$'Lat.50%')))
  proposal.z <- mvnorm(S = diag(c(0.0025, 0.0025)),
                       n = nlocation(cbind(zsum$'Lon.50%', zsum$'Lat.50%')))
  
  cat("Running final model for ", BirdId[i], ".........\n")
  
  # Note the increase in number of interations - this takes a bit longer to run
  FIT[[i]] <- estelleMetropolis(
    model = model,
    proposal.x = mvnorm(chainCov(fit$x), s = 0.1),
    proposal.z = mvnorm(chainCov(fit$z), s = 0.1),
    x0 = cbind(xsum$'Lon.50%', xsum$'Lat.50%'),
    z0 = cbind(zsum$'Lon.50%', zsum$'Lat.50%'),
    iters = niters,
    # This value sets the number of iterations to run
    thin = nthin,
    chains = 3
  )
  cat("\n\n\n")
}

# saveRDS(FIT, file = "D:/OSFL_Geolocator/osfl_fit_Mar_2021_AbundancePrior_AncZenith.rds")

Export daily data for movebank

FIT <-
  readRDS("G:/OSFL/osfl_fit_Feb_2021_AbundancePrior_AncZenith.rds")

# BirdId <- c(BirdId, "K288_yr2")

movebank <- vector('list', (nBirds + 1))
for (i in 1:(nBirds + 1)) {
  movebank[[i]] <-
    SGAT2Movebank(FIT[[i]]$z, time = FIT[[i]]$model$time)
  movebank[[i]]$TagID <- BirdId[i]
}

# write.csv(do.call('rbind', movebank),row.names = FALSE, "F:/OSFL/movebank_export_Mar24_2021.csv")

saveRDS(FIT, paste0("OSFL_",Sys.Date()))